1 Descripción del ejercicio

El presente ejercicio se engloba dentro de la materia de Sistemas de Información Geográfica y Ecología Espacial: Aplicaciones del Máster Geoforest de la Universidad de Córdoba.

El trabajo consiste en aplicar las herramientas aprendidas en clase, con lo que los productos finales esperados son:

1.1 Datos de partida

Para la realización del ejercicio partimos de los siguientes datos:

  • limite_yunquera.shp: límite de la zona de estudio (EPSG: 25830).

  • vegetacion_yunquera.shp: mapa de vegetación de la zona de estudio.

  • inventario_pinsapo.shp: puntos de inventario de parcelas de pinsapo (Abies pinsapo).

  • PNOA_MDT05_ETRS89_HU30_1051_LID.asc: modelo digital del terreno con tamaño de píxel de 5x5 metros.

Estos datos se encuentran en la carpeta Cartografia/CartografiaOriginal.

El resto de datos de partida se obtendrán a lo largo de la resolución del ejercicio.

1.2 Zona de estudio

Antes de comenzar a resolver el ejercicio conviene hacer una breve descripción de la zona de estudio.

Yunquera es un municipio perteneciente a la provincia de Málaga (Andalucía). En este monte se encuentra una parte importante del pinsapo. En el siguiente mapa se pueden explorar las diferentes formaciones vegetales presentes en el área de estudio (hacer click en las parcelas para explorar).


2 Resolución del ejercicio

Se ha dividido el ejercicio en 4 bloques de acuerdo a los objetivos del mismo.

2.1 Flujo de trabajo

2.2 Mapa de biodiversidad

En esta sección se construirá un mapa donde se representará el Índice de Shanon:

\[ \tag1 H' = - \sum^s_{i=1} p_i \times log_2 p_i \]


donde:

  • i: cada una de las distintas especies

  • s: total de especies distintas

  • \(p_i\): abundancia relativa de cada especie en la comunidad (N individuos de la especie i / N total de individuos).

2.2.1 Descarga de ocurrencias

En primer lugar debemos descargar las ocurrencias de GBIF. En este caso se ha creado un rectángulo en la plataforma web y se han extraído todas las ocurrencias que se encontraban dentro de este rectángulo (correspondiente a la zona de estudio).

A continuación se importan las ocurrencias en R con el nombre de puntos. El siguiente paso es convertir estos puntos en un objeto espacial. Para ello se utiliza la función sf::st_as_sf(). La salida de esta operación (puntos_sf) se transforma al sistema de coordenadas en el que se está trabajando (EPSG: 25830) llamando al objeto ocurrencias. En la Fig. 1 comprobamos que hemos elegido correctamente el rectángulo, y que nuestros puntos se extienden más allá de la zona de estudio asegurando la inclusión de todas las ocurrencias.

puntos <- read.csv(file = here('Cartografia/Tablas/ocurrences.csv'),
                   sep = "\t",
                   quote = "") |> 
  dplyr::select(scientificName,species,decimalLatitude,decimalLongitude) |> 
  drop_na()

puntos_sf <- st_as_sf(puntos,
                      coords = c('decimalLongitude','decimalLatitude'),
                      crs = 4326)

ocurrencias <- st_transform(puntos_sf,
                            crs = 25830)
<center>Fig. 1. Localización de los puntos descargados de GBIF con respecto al área de estudio</center>

Fig. 1. Localización de los puntos descargados de GBIF con respecto al área de estudio


El siguiente paso será delimitar los puntos escogidos al rectánculo mínimo envolvente (bounding box) del área de estudio mediante la función sf::st_crop(). La razón de escoger todos los puntos incluidos dentro del rectángulo y no solamente los que están dentro del área de estudio es que como veramos en la siguiente sección trabajaremos con comunidades ecológicas de forma cuadrangular.

ocurrencias <- st_crop(ocurrencias, yunquera)

Podemos ver el resultado en la Fig. 2, que ascienden a un total de 16906 ocurrencias correspondientes a 1465 especies diferentes.

<center>Fig. 2. Ocurrencias de GBIF dentro del monte de Yunquera</center>

Fig. 2. Ocurrencias de GBIF dentro del monte de Yunquera


2.2.2 Creación de comunidades ecológicas

Las comunidades ecológicas se asumirán que son entidades cuadrangulares de 250 m de lado. Para ello creamos primero la base sobre la que se crearán estas comunidades ecológicas, que será el retángulo mínimo envolvente de la zona de estudio (bounding box), el cuál se construye con la función sf::st_bbox(). A continuación se crea una malla (sf::st_make_grid()) con tamaño de celda de 250x250m (cellsize), y se les asigna un GID único a cada una. Finalmente se crea el objeto ocurrencias_id que no es más que la misma entidad de ocurrencias que teníamos con un campo nuevo correspondiente al GID de la comunidad ecológica a la que pertenecen. La última función crea un data frame llamado bio que será con el que calcularemos el índice de Shannon.

# Crear bounding box y convertir a espacial
BBox <- st_bbox(yunquera)
BBox <- st_as_sfc(BBox)

# Crear malla de 250x250m (utilizar old pipeline para nchar)
grid250  <- BBox %>%
  st_make_grid(square = TRUE, cellsize = c(250, 250)) %>%
  cbind(data.frame(ID = sprintf(paste("GID%0",
                                      nchar(length(.)),
                                      "d",
                                      sep=""), 
                                1:length(.)))) %>%
  st_sf()

# Añadir ID de comunidad ecologica a cada ocurrencia
ocurrencias_id <- ocurrencias |> 
  st_join(grid250, join = st_intersects ) |> 
  mutate(id_250 = factor(ID),
         scientific = scientificName) |> 
  dplyr::select(-ID,-scientificName)

# Crear data frame 
bio <- data.frame(id_250 = ocurrencias_id$id_250,
                  scientific = ocurrencias_id$scientific) |> 
  drop_na()
Tabla 1. Encabezado de la tabla de datos de ocurrencias en las distintas cuadrículas
id_250 scientific
GID057 Bupleurum spinosum Gouan & Ill.Observ.Bot.
GID107 Phoenicurus ochruros (S.G.Gmelin, 1774)
GID026 Phlomis crinita Cav.
GID026 Podarcis vaucheri (Boulenger, 1905)
GID026 Quercus faginea Lam.
GID105 Abies pinsapo Boiss.


En la Fig. 3 tenemos una representación gráfica de como quedaría esta malla cuadrangular con respecto a las ocurrencias y a la zona de estudio, y en la Fig. 4 vemos la densidad de puntos en la zona.

<center>Fig. 3. Malla cuadrangular de 250 x 250 m en el monte de Yunquera. Los puntos son todas las ocurrencias consideradas en los análisis posteriores</center>

Fig. 3. Malla cuadrangular de 250 x 250 m en el monte de Yunquera. Los puntos son todas las ocurrencias consideradas en los análisis posteriores


<center>Fig. 4. Mapa de densidad de ocurrencias en el monte de Yunquera</center>

Fig. 4. Mapa de densidad de ocurrencias en el monte de Yunquera


2.2.3 Cálculo del índice de Shannon

El último paso de esta sección será calcular el Índice de Shannon (H). Esto lo haremos en cuatro pasos:

  1. Calcular el número de especies por cuadrícula. En el siguiente cuadro de código se agrupa por cuadrícula y nombre (group_by), se genera una columna con el número de especies (count()), y se le cambia el nombre a la columna. En la Tabla 2 vemos el encabezado del resultado.
T_num_ind_sp_cuad <- bio |> 
  group_by(id_250, scientific) |> 
  count() |> 
  rename(num_ind_sp_cuad = n)
Tabla 2. Encabezado de la tabla del número de especies por cuarícula
id_250 scientific num_ind_sp_cuad
GID001 Erinacea anthyllis Link 3
GID001 Quercus faginea Lam. 6
GID001 Thymus mastichina subsp. mastichina 3
GID002 Erinacea anthyllis Link 6
GID002 Quercus faginea Lam. 12
GID002 Thymus mastichina subsp. mastichina 6


  1. Calcular el número de ocurrencias por cuadrícula. Realizamos un procedimiento similar, pero en este caso solamente agrupamos por cuadrícula. Vemos el encabezado en la Tabla 3.
T_num_ind_cuad <- bio |> 
  group_by(id_250) |> 
  count() |> 
  rename(num_ind_cuad = n)
Tabla 3. Encabezado de la tabla del número de ocurrencias por cuadrícula
id_250 num_ind_cuad
GID001 12
GID002 24
GID003 53
GID004 37
GID005 17
GID006 25


  1. El siguiente paso calcula el Índice de Shannon (T_Shannon). Para ello se unen las tablas haciendo un left_join() que incluye todas las columnas de la Tabla 1 de forma que cada cuadrícula tendrá asociada sus especies, el número de ocurrencias de cada una, y el número total de ocurrencias de la cuadrícula de esa especie. A continuación se crean los coeficientes \(p_i\) y \(\log p_i\) que se utilizarán para el cálculo del índice. En la Tabla 4 vemos la estructura de los datos hasta este punto.
    Finalmente, se calcula el índice de Shannon, se seleccionan las columnas necesarias y se eliminan los duplicados.
T_Shannon <- T_num_ind_sp_cuad |> 
  left_join(T_num_ind_cuad,
            by = 'id_250') |> 
  mutate(pi = num_ind_sp_cuad/num_ind_cuad,
         logpi = log2(pi)) |> 
  group_by(id_250) |> 
  mutate(shannon = sum(logpi * pi)*-1) |> 
  dplyr::select(id_250, shannon) |> 
  distinct()
Tabla 4. Coeficientes para el cálculo del Índice de Shannon
id_250 scientific num_ind_sp_cuad num_ind_cuad pi logpi
GID001 Erinacea anthyllis Link 3 12 0.25 -2
GID001 Quercus faginea Lam. 6 12 0.50 -1
GID001 Thymus mastichina subsp. mastichina 3 12 0.25 -2
GID002 Erinacea anthyllis Link 6 24 0.25 -2
GID002 Quercus faginea Lam. 12 24 0.50 -1
GID002 Thymus mastichina subsp. mastichina 6 24 0.25 -2


  1. El último paso consiste en dar el valor de Índice de Shannon a cada una de las cuadrículas y reemplazar donde no hay ninguna ocurrencia por un 0. Vemos en el mapa la representación final
# Valor H a cada cuadricula
SF_shannon <- grid250 |> 
  left_join(y = T_Shannon,
            by = c("ID" = "id_250"))

# Sustituir NA por 0
SF_shannon$shannon <- replace_na(SF_shannon$shannon, 0)

# Corte para visualización
SF_shannon_inter <- st_intersection(SF_shannon, yunquera)


2.3 Segmentación

La segmentación consiste en dividir el territorio en unidades homogéneas en cuanto a unas determinadas características de las variables de entrada. El procedimiento consiste de los siguientes pasos:

  1. Elegir las variables de entrada

  2. Convertir a formato raster y homogeneizar las escalas

  3. Construir un raster virtual

  4. Segmentación del territorio

2.3.1 Variables de entrada

En esta sección se realizarán los pasos 1 y 2 previos. Para ello, construiremos un total de 4 variables (ver cuadro siguiente).

Variables de entrada

Densidad de pinsapo: información sobre la vegetación

Orientaciones: información sobre el estrés abiótico derivado de la situación de solana y umbría

Distancia a canales: información sobre zonas de posible compensación hídrica

Índice de Shannon: información sobre la biodiversidad


Densidad de pinsapo

En primer lugar cargamos la capa, la cual se encuentra en CRS 25830, y la visualizamos con los límites del monte para comprobar que todo es correcto.

SF_pinsapo <- read_sf(here("Cartografia/CartografiaOriginal/inventario_pinsapo.shp"),
                      crs = 25830)


El siguiente paso será crear un raster de densidad. En ese sentido, se utilizará la interpolación inversa a la distancia (IDW).

  1. Primero se transforman los puntos a WGS 1984 ya que las funciones siguientes trabajan este sistema de coordenadas.
SF_pinsapo84 <- st_transform(SF_pinsapo, crs = 4326)
  1. A continuación obtenemos el rectángulo mínimo envolvente(BBox84). Se da una toleracia de 0.005º dado que al reproyectar la imagen se pierde información.
BBox84 <- st_bbox(SF_pinsapo84)


BBox84[c(1,3)] <- BBox84[c(1,3)] + c(-0.005,0)
BBox84[c(2,4)] <- BBox84[c(2,4)] + c(-0.005,0)
  1. Creamos un patrón de puntos con la función spatstat.geom::ppp() utilizando el vector de coordendas X e Y, marks es la variable que utilizaremos para interpolar (Nparc), y en un objeto owin() ponemos las coordendas X e Y del rectángulo mínimo envolvente, que será la extensión de salida.
ppp_pinsapo <- ppp(st_coordinates(SF_pinsapo84)[,1],
                   st_coordinates(SF_pinsapo84)[,2],
                   marks = SF_pinsapo$Nparc,
                   window = owin(BBox84[c(1,3)],
                                 BBox84[c(2,4)]))
  1. Con la función spatstat::idw() realizamos la interpolación y se transforma a raster, que en este caso a diferencia de QGIS no tenemos opción de escoger el número de vecinos de cuáles hacer la interpolación (no obstante, la diferencia del raster de salida es mínima).
R_pinsapo <- ppp_pinsapo |> 
  idw(power = 2, at = 'pixels') |> 
  raster()
  1. Se asigna la proyección y se reproyecta al sistema 25830 cambiando la resolución a 10 metros.
# Asignar proyección
crs(R_pinsapo) <- crs(SF_pinsapo84)

# Transformar a 25830 y asignar resolución espacial
R_pinsapo <- projectRaster(R_pinsapo,
                           crs = 25830,
                           res = 10)
  1. Finalmente se hace una reclasificación de los valores del ráster y se recorta al área de estudio.
# Matriz de entrada para reclasificación
classMatrix <- matrix(c(0,50,0,
                        50,200,100,
                        200,400,300,
                        400,600,500,
                        600,800,700,
                        800,Inf,900),
                      byrow = T,
                      ncol = 3)

# Reclasificación y recorte al área de estudio
R_pinsapo <- R_pinsapo |> 
  reclassify(rcl = classMatrix) |> 
  mask(yunquera)

En el siguiente webmap podemos ver el resultado de los análisis anteriores.


Orientaciones

La segunda de las capas se corresponde con la de orientaciones. En este caso cargamos el modelo digital del terreno (mdt), calculamos las orientaciones en grados con la función raster::terrain(). El siguiente paso será hacer coincidir espacialmente ambos raster con la misma resolución (resample()). Finalmente lo reclasificaremos teniendo en cuenta que valores bajos y valores altos significan lo mismo (Norte) y lo trataremos como un raster con valores discretos con la función ratify().

# Cargar capa y cortar a la extensión de Yunquera
mdt <- raster(here("Cartografia/CartografiaOriginal/PNOA_MDT05_ETRS89_HU30_1051_LID.tif"),
              crs = "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs") 

# Cálculo orientaciones
orientaciones <- terrain(mdt,
                         opt = 'aspect',
                         unit = 'degrees') |> 
  resample(R_pinsapo)

# Matriz de clasificación
mat <- matrix(c(-Inf, 22.5,1,
                22.5,67.5,2,
                67.5,112.5,3,
                112.5,157.5,4,
                157.5,202.5,5,
                202.5,247.5,6,
                247.5,292.5,7,
                292.5,337.5,8,
                337.5,Inf,1),
              ncol = 3,
              byrow = TRUE)

# Reclasificación y recorte al área de estudio
R_orientaciones <- orientaciones |> 
  reclassify(rcl = mat) |> 
  mask(yunquera) |> 
  ratify()


Distancia a canales

El ráster de distancia a canales se ha calculado previamente siguiendo los pasos del guión en QGIS debido a la complicidad de realizar este análisis en R. La capa obtenida se muestra en el siguiente webmap, la cual se encuentra en el sistema de coordenadas correspondiente y con resample() nos aseguramos de que coincida espacialmente con el resto de capas.

R_dist_canales <- raster(here("Cartografia/CartografiaOriginal/Raster_distancia_cauces.tif")) |> 
  resample(R_pinsapo) |> 
  mask(yunquera)


Indice de Shannon

El indice de Shannon se ha calculado en la sección 2.2.3. En este caso nos quedaría rasterizar el resultado y remuestrearlo para que coincida espacialmente con el resto de capas.

R_shannon <- SF_shannon |>
  dplyr::select(shannon) |> 
  st_rasterize() |> 
  as("Raster") |> 
  resample(R_pinsapo) |> 
  mask(yunquera)

Guardar las capas

Finalmente se exportan las capas a la carpeta CapasDefRaster.

# writeRaster(R_pinsapo,
#             filename = here("Cartografia/CapasDefRaster/R_pinsapo.tif"),
#             overwrite = T)
# 
# writeRaster(R_orientaciones,
#             filename = here("Cartografia/CapasDefRaster/R_orientaciones.tif"),
#             overwrite = T)
# 
# writeRaster(R_dist_canales,
#             filename = here("Cartografia/CapasDefRaster/R_dist_canales.tif"),
#             overwrite = T)
# 
# writeRaster(R_shannon,
#             filename = here("Cartografia/CapasDefRaster/R_shannon.tif"),
#             overwrite = T)


2.3.2 Raster virtual

El ráster virtual se ha creado en QGIS mediante la herramienta Raster -> Micelanea -> Crear raster virtual.

El resultado se muestra en la Fig. 5.

<center>Fig. 5. Creación del Raster virtual en QGIS a partir de las cuatro capas de entrada</center>

Fig. 5. Creación del Raster virtual en QGIS a partir de las cuatro capas de entrada


2.3.3 Segmentación del territorio

El último paso de esta sección consiste en realizar una segmentación del monte de Yunquera usando el ráster virtual creado en la sección anterior. Para ello, se utilizará la herramienta Orfeo Toolbox en la interfaz de QGIS.

Para llevar esta tarea a cabo, con Orfeo correctamente instalado, se debe ir a Caja de herramientas de Procesos -> OTB -> Segmentation -> Segmentation. Una vez ahí, se seleccionaron las opciones de la Fig. 6. La opción “Minimum region size” se probó con varios valores. En el webmap que se verá posteriormente se analizarán las diferencias e influencia de las diferentes capas de entrada en los resultados.

<center>Fig. 6. Parámetros utilizados en la segmentación</center>

Fig. 6. Parámetros utilizados en la segmentación


Para finalizar la segunda parte de la práctica, se cargan las tres capas creadas con Orfeo en QGIS.

segmentation1500 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion1500.shp"))

segmentation2000 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion2000.shp"))

segmentation2500 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion2500.shp"))

En el webmap siguente se puede jugar con las capas para ver el efecto de las distintas segmentaciones en cada una de las capas.


Nota: los bordes de la segmentación no coinciden totalmente con el borde de los ráster. Esto se debe a que las capas se encuentran en crs 25830, pero el webmap se representa en WGS 1984.

Los resultados se discuten en la última sección. No obstante, Los resultados se harán en base a la segmentación escogida que se discute en las siguientes líneas.

Para analizar los resultados podemos activar en primer lugar Segmentacion (2000), y a continuación Segmentacion (2500) (el orden es importante). Tras hacer esto, podemos activar las variables una a una empezando por Densidad pinsapo, Indice de Shannon, Orientaciones y Distancia a canales (m). La segmentación de 2500 píxeles como unidad mínima se tratará como segmentación inicial, y las unidades segmentadas como cantones.

Lo primero que podemos ver es que si segmentamos con unidad mínima de 2000, se crean 7 cantones a mayores. La segmentación inicial parece separar muy bien áreas con alta densidad de pinsapo de áreas con baja densidad. Al aumentar el número de cantones no se ve una clara separación de nuevas áreas con densidad de pinsapo diferentes. Si activamos el índice de Shannon vemos que no existe una clara relación entre los cantones y la diver


2.4 Webmap

El tercer objetivo era crear un webmap y publicarlo. En este caso se ha publicado en GitHub el webmap previo. Se puede ver en el siguiente enlace.

2.5 Cuestiones finales

  • Consideras que el resultado es consistente con la realidad que conoces del terreno? Justifica tu respuesta.